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ABSTRACT 

Intensity mapping, which images a single spectral line from unresolved galaxies across cosmological volumes, 
is a promising technique for probing the early universe. Here we present predictions for the intensity map 
and power spectrum of the CO(l-O) line from galaxies at z ~ 2.4-2.8, based on a parameterized model for 
the galaxy-halo connection, and demonstrate the extent to which properties of high-redshift galaxies can be 
directly inferred from such observations. We find that our fiducial prediction should be detectable by a realistic 
experiment. Motivated by significant modeling uncertainties, we demonstrate the effect on the power spectrum 
of varying each parameter in our model. Using simulated observations, we infer constraints on our model 
parameter space with an MCMC procedure, and show corresponding constraints on the L\v-Lco relation and 
the CO luminosity function. These constraints would be complementary to current high-redshift galaxy obser¬ 
vations, which can detect the brightest galaxies but not complete samples from the faint end of the luminosity 
function. By probing these populations in aggregate, CO intensity mapping could be a valuable tool for probing 
molecular gas and its relation to star formation in high-redshift galaxies. 

Subject headings: galaxies: high-redshift — galaxies: evolution — ISM: molecules 


1. INTRODUCTION 

Current progress in cosmology and galaxy formation is 
strongly informed by galaxy observations at increasingly high 
redshifts, reaching into the epoch of reionization. Modern ob¬ 
servatories such as ALMA may have even begun to probe the 
interstellar medium (ISM) of “typical” high-redshift galaxies 
(e.g. Riechers et al. 2014; Watson et al. 2015). Despite this, 
it remains difficult and expensive to observe populations of 
such galaxies in statistically large and complete samples, due 
to their faintness. 

However, a full understanding of both cosmology and 
galaxy formation depends on such samples. To inform galaxy 
formation models, individual galaxies need to be placed in the 
context of their broader populations, for without such context 
it is unclear how representative individually observed galaxies 
may be. For cosmology, precise clustering measurements that 
trace large-scale stmcture require sufficient number densities 
of spectroscopically observed galaxies (Blake et al. 2011; An¬ 
derson et al. 2012). 

Intensity mapping is a technique that can potentially ad¬ 
dress this challenge in regimes that are inaccessible by typical 
galaxy surveys. The method aims to map a single spectral 
line across cosmologically large volumes. The angular reso¬ 
lution used to do this images cumulative emission from multi¬ 
ple galaxies, instead of resolving single galaxies. In this way, 
intensity mapping aggregates flux from galaxies that would 
individually be below the detection limit, while still resolv¬ 
ing large-scale stmcture by measuring intensity fluctuations 
over cosmological distances. As a consequence, this tech¬ 
nique will necessarily describe galaxies in a statistical sense 
through measures like the spatial power spectrum. 

In this paper, we look at intensity mapping of carbon 
monoxide (CO), a common tracer of molecular gas and star 
formation in nearby galaxies. After H 2 , CO is the most 


abundant molecular species, tracing the metal-enriched, rela¬ 
tively dense (> 10 2 cm -3 ), cool to warm molecular ISM phase 
where stars form efficiently. This physically motivates the 
empirical conversion between CO and H 2 (for a review, see 
Bolatto et al. 2013), and by extension star formation. The CO 
lines themselves arise from a “ladder” of rotational transitions 
starting at millimeter wavelengths. For this study, we focus on 
the ground-state CO(l-O) transition at 115.27 GHz (2.6 mm). 

Other lines besides CO have also been proposed as can¬ 
didate lines for intensity mapping. The Hi 21 cm line, in 
particular, has been extensively studied for the purposes of 
tracing large-scale structure out to z ~ 2.5 as well as imaging 
hydrogen reionization at z > 6 (Chang et al. 2008; Morales & 
Wyithe 2010; Bandura et al. 2014). Additional potential tar¬ 
gets include [Cll] at 158 pm (Gong et al. 2012; Silva et al. 
2015; Uzgil et al. 2014), Lya at 1216 A (Pullen et al. 2014), 
and various other fine structure lines (Visbal & Loeb 2010), 
each with its own advantages and challenges. 

Observing CO—as well as Hi, [Cll], and other lines—is 
necessary for a census of all phases of the ISM in galaxies. 
However, a practical advantage of CO is that it simultaneously 
emits at multiple frequencies: because of the small excitation 
energies of the lowest energy levels (5.5 and 16.6 K for the 
/ = 1 and 2 levels, respectively), a galaxy is likely to have 
comparable emission in CO(l-O) and CO(2-l). Observed sig¬ 
nals in carefully chosen frequency bands can potentially be 
cross-correlated to determine the contribution of the CO-only 
signal. 

Our focus is the “epoch of galaxy assembly,” at redshifts of 
roughly z ~ 2-3. This is a particularly interesting epoch for 
galaxy formation, because it is near the peak of cosmic star 
formation. Current and near-future galaxy surveys are also ex¬ 
pected to reach into these redshifts, allowing the opportunity 
to cross-correlate the CO signal with galaxy surveys (Pullen 
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Fig. 1.— Input and output of our modeling process, i.e. initial dark matter halos and final CO intensity map (details in §2.2). These plots illustrate one 
realization of the pathfinder experiment’s survey volume (§2.4 and Table 2), while the full experiment’s survey area is 2.5 times larger. Top: Halos in the 
3D volume, rendered to scale in comoving distance. Along the line-of-sight direction, we label the equivalent cosmological redshifts and redshifted CO(l-O) 
frequencies. Middle: 2D projections of halo positions. The left image shows the “front” view of halos that would fall into the highest 40 MHz frequency channel, 
or lowest redshift slice. The pathfinder beam size is shown for scale. The right image shows the “side” view of halos to a depth of 6 arcmin, or one beam 
width. Bottom: CO intensity map produced by our fiducial model. The slice volumes are the same as above, albeit with comoving depth converted to observed 
frequency. The same large-scale structure is readily apparent in both images, even with the lower resolution of the intensity map. The analysis in this paper relies 
on the power spectrum of this map (see Fig. 3). 


et al. 2013). Our current understanding of star formation and 
gas content in this epoch is incomplete, and largely limited 
to the bright end of the relevant populations. In the longer 
term, observations at these redshifts could serve as a stepping 
stone for future CO observations that reach into the epoch of 
reionization (Carilli 2011; Gong et al. 2011; Lidz et al. 2011). 

Previous predictions for the intensity of the CO signal vary 
by more than an order of magnitude (Breysse et al. 2014, at 
z ~ 3). The wide range simply reflects the current scarcity 
of data for typical high-redshift galaxies. It is possible to di¬ 
rectly simulate these galaxies, but such simulations are expen¬ 
sive and still are quite uncertain. These uncertainties suggest 


a need for alternative probes of high-redshift galaxy popula¬ 
tions, especially over numbers and/or volumes currently inac¬ 
cessible to traditional surveys. 

Given the modeling uncertainties, predictions of the ex¬ 
pected signal will only go so far, at least until a measurement 
is attempted. Here we also ask, what could we learn from in¬ 
tensity mapping if a measurement is made? More precisely: 
given hypothetical but tractable intensity mapping observa¬ 
tions, what can we infer about the properties and distribution 
of the underlying galaxy population? To our knowledge, this 
question has not yet been directly addressed in the literature. 
Here we put these questions in the context of CO surveys that 
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Fig. 2.— Relation between mean CO luminosity and halo mass at z = 2.4, 
at the low end of the redshift range in this study. The black solid line shows 
the relation used in this work. For comparison, Lco(M/,) from previous stud¬ 
ies (Righi et al. 2008; Visbal & Loeb 2010; Pullen et al. 2013) have also been 
plotted 2 . Solid color lines indicate the exact linear Lco(M/,) relations used 
in those studies. Dotted color lines indicate approximate linear scalings for 
models without a direct Lco(M/,) relation (as in Breysse et al. 2014). In our 
model, we imposed a minimum CO-luminous halo mass of 10 10 M@, exclud¬ 
ing halos in the shaded gray range (see Appendix A for some justification). 

a To compare models consistently, the Lco(M /,) relations plotted here 
have absorbed a “duty-cycle ” factor, calculated in those papers as /d ut „ ~ 
10 s yr/f a ge(z). Carilli (2011), Gong et al. (2011), and Lidz et al. (2011) 
are omitted from this plot because they focus only on reionization redshifts 
(z>6). 

are now being built or planned for the near future. 

The structure of this paper is as follows. In §2, we sum¬ 
marize our method for (1) modeling galaxy CO luminosities 
and generating CO intensity maps and power spectra and (2) 
inferring constraints on this model given a hypothetical obser¬ 
vation. In §3, we show how the variation of different model 
parameters affects our fiducial prediction. We then present 
MCMC-inferred parameter constraints from mock intensity 
mapping observations, as well as the inferred Lir-Lco relation 
and CO luminosity function. In §4, we discuss our results, 
the uncertainties in our modeling, and the ability of intensity 
mapping to constrain the properties of galaxy populations. Fi¬ 
nally, in §5, we resummarize our main results and conclude. 

2. METHODS 

2.1. Context: Previous Studies 

Several previous studies have explored CO intensity map¬ 
ping and modeled the expected signal. To date, they include: 
Righi et al. (2008), Visbal & Loeb (2010), Gong et al. (2011), 
Carilli (2011), Lidz et al. (2011), Pullen et al. (2013). 

Since the models in those studies inform our own, the 
reader may find a brief summary useful. Table 3, located to¬ 
ward the end of this paper, compares the most relevant aspects 
of each model. The common thread across nearly all of them 
(except Gong et al. 2011) is an assumed connection between 
star-formation rate (SFR) and CO luminosity. In those papers, 
the procedure is to first calculate SFRs and then convert to CO 
luminosity, though the details vary. 

Around the redshifts considered in this study, Breysse et al. 
(2014) found that the mean CO intensity predicted by some of 
these models differs at the extremes by more than an order of 
magnitude, which simply reflects the significant uncertainty 


in these predictions. 

2.2. Part 1: Modeling CO luminosity 

In our modeling, we begin with dark matter halos (tracing 
underlying cosmological structure) and ultimately generate 
simulated three-dimensional CO intensity maps. 

In order, the components of this calculation are: 

1. Dark matter halos (§2.2.1) 

2. Star-formation rates (§2.2.2) 

3. Infrared luminosities (§2.2.3) 

4. CO luminosities (§2.2.4) 

5. Intensity maps and power spectra (§2.2.5) 

Figure 1 illustrates the input (dark matter halos) and output 
(CO intensity map) of this process. It also provides some vi¬ 
sual intuition for the dimensions of the survey volume, as well 
as the cosmological structure contained within. 

We discuss each of the individual steps in more detail be¬ 
low. Where parameters are introduced, we state their fiducial 
values, as well as the priors adopted for the MCMC analysis 
(§2.3). 

2.2.1. Dark Matter Halos 

We begin with dark matter halos in a “lightcone" volume. 
To obtain this, we use the results of a cosmological /V-body 
dark matter simulation (the c400-2048 box 1 ). The simula¬ 
tion was run with L-Gadget (based on Gadget- 2, Springel 
et al. 2001; Springel 2005). The box has 2048 3 particles and 
a side length of 400 Mpc hr 1 , resulting in a particle mass of 
5.9 x 10 8 M Q /T I . The softening length used is 5.5 kpc hr 1 , 
constant in comoving length. The initial conditions are gen¬ 
erated by 2LPTIC 2 (Crocce et al. 2006) at z = 99, with the 
power spectrum generated by Camb 3 . 

The ACDM cosmological parameters of the simulation are 
n,„ = 0.286, n A = 0.714, il h = 0.047, h = 0.7, erg = 0.82, and 
n s = 0.96. We assume these values throughout this paper. 

Dark matter halos were identified at each simulation snap¬ 
shot using the Rockstar halo finder (Behroozi et al. 2013c). 
In this work, we treated subhalos and central halos identically, 
but excluding subhalos from our calculations does not signif¬ 
icantly affect our results. 

The simulation output consisted of 100 snapshots, from z = 
12.33 to z = 0 inclusive, with equal logarithmic spacing in 
(1+z) -1 , though only the halos from six of these snapshots 
were ultimately used. Lightcones were generated from these 
snapshots by choosing an arbitrary z = 0 observer origin and 
direction, then selecting all halos along the line of sight within 
the desired survey area and redshift range (see §2.4) from the 
appropriate redshift snapshots. 

2.2.2. Halos and Star Formation 

To get the SFR of a dark matter halo, we use the results 
of Behroozi et al. (2013a) and Behroozi et al. (2013b), which 
empirically quantified the stellar mass history of dark matter 

1 Provided by Matthew Becker (M. Becker et al. 2015, in preparation) 

2 http : //cosmo . nyu . edu/roman/ 2 LPT/ 

3 http : //camb . info/ 
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TABLE 1 

Model parameters 


Label 

Link 

Description 

Fiducial 

Prior 

Details 

O'SFR 

Halos —> SFR 

Log-scatter in SFR 

0.3 

A/(0.3,0.1) and ctsfr > 0 

§2.2.2 

log <5 mf 

SFR — > Lir 

SFR-Ljr scaling 

0.0 

Af (0.0,0.3) 

§2.2.3, Eq. 1 

a 

Lir — > Leo 

Lir—L eo log-slope 

+1.37 

Af(1.17,0.37) 

§2.2.4, Eq. 2 

p 

L\r — > L C o 

Lir-Lco log-intercept 

-1.74 

Af(0.21,3.74) 

§2.2.4, Eq. 2 

a Lco 

Lir -a L C o 

Log-scatter in Leo 

0.3 

A/(0.3,0.1) and cr Lm >0 

§2.2.4 


Note. — A7//. ct) denotes a normal distribution with mean // and variance ct 2 . 


TABLE 2 

Instrument and Survey Parameters 


Type 

Description 

Pathfinder 

Full 

Instrument 

System temperature (L^s) 

40 K 

35 K 


Dual-polarization feeds (Nf ee ds) a 

19 

500 


Beam width (#fwhm) 

6' 

3' 


Frequency band (A v) 

30-34 GHz 

30-34 GHz 


Frequency channels (<5^) 

40 MHz 

10 MHz 

Survey 

Survey area / patch (f2 S urv) 

2.5 deg 2 

6.25 deg 2 


On-sky time / patch (r 0 b s ) b 

1500 hr 

2250 hr 

Derived 

CO(l-O) redshift range 

OO 

ri 

1 

<N 

2.4-2.8 


Instrument sensitivity 

1026 fiK s 1 / 2 

783 iiK s 1 / 2 


Final map sensitivity c 

41.5 ix K MHz 1 / 2 

18.3 //.K MHz 1 / 2 


Total P(k ) detection significance (SNR to t) d 

7.89 ct 

144ct 


“Dual polarization effectively doubles the number of feeds, so in our calculations we model (Vf ee ds = 19 1500) as 38 (1000) feeds. 

b We assume that observing time is divided equally over four patches with ~35% observing efficiency. Thus, 750 hr per patch equals a total of ~1 year of 
operation. 

“This sensitivity measure is proportional to the 3D noise power spectral density (Eq. C2). It is independent of spectral resolution, which is important here since 
we are mainly interested in measuring 3D structure. To get rms noise per channel instead, divide the value by 
d See §3. 1 for additional discussion. 


halos back to z = 8. Briefly, that study constrained a parame¬ 
terized stellar mass-halo mass relation by applying it to sim¬ 
ulated halo merger trees, accounting for systematic errors and 
biases, and comparing derived stellar mass functions, cosmic 
SFRs, and specific SFRs with a comprehensive compilation 
of observational data. The reader is referred to the original 
papers for a more detailed discussion of that study. 

For our purposes, we are primarily interested in SFR(M,z), 
their derived results for the average halo SFR for a given halo 
mass and redshift. To simplify this calculation, we interpolate 
their tabulated data 4 for SFR(M,z). 

The results of Behroozi et al. (2013a) also provide ±lcr 
posterior constraints on SFR. Note that this is distinct from 
halo-to-halo scatter, which we address later. However, be¬ 
cause these bounds span a fairly consistent log-space inter¬ 
val (~ ±0.15 dex) over the relevant masses and redshifts, any 
variation within them is effectively a rescaling factor, which 
we absorb into the parameter (Smf (see below, §2.2.3, Eq. 1). 

We also add halo-to-halo log-normal scatter, parameterized 
as ctsfr, since the results of Behroozi et al. (2013a) constrain 
only average halo SFRs. This single-parameter scatter is a 
simple way of capturing the variation in SFR for a given halo 
mass. There is evidence that “normal” star-forming galaxies 
exhibit a strong correlation with their stellar mass, with a scat¬ 
ter of ~ 0.2-0.4 dex, while starbursts may exist as a separate 
population with unusually high SFR (e.g. Speagle et al. 2014; 
Salmon et al. 2015). We assume the scatter in SFR given stel¬ 
lar mass to be reasonably similar to the scatter in SFR given 
halo mass. With this in mind, we choose a fiducial value of 
csfr = 0.3 and a prior of ctsfr = 0.3 ±0.1. To avoid unphysical 

4 Available at the time of writing at http : //www . peterbehroozi . 
com/data . html 


negative scatter, we also require ctsfr > 0. 

2.2.3. Star Formation and Infrared Luminosity 

Empirically connecting SFR to Leo requires at least one in¬ 
termediate, directly observable quantity, since SFRs are not 
directly measured. Instead, they are inferred from photomet¬ 
ric or spectral tracers with various underlying assumptions. 
We will take this intermediate tracer to be the total infrared 
luminosity Lir (see Carilli 2011; Lidz et al. 2011; Pullen et al. 
2013), conventionally the integrated 8-1000 fim luminosity. 
Thus, the model is calibrated on empirical correlations be¬ 
tween SFR, Lir, and Leo, resting on physical assumptions 
about star formation, dust, and molecular gas. Some recent 
work has focused on quantifying these relations at higher 
redshifts (e.g. Bethermin et al. 2015, for z < 4), but signif¬ 
icant uncertainties remain, especially in fainter galaxies and 
at higher redshifts, and part of our aim is to analyze how the 
signal may vary within those uncertainties. 

We assume a correlation between SFR and Lir (Kennicutt 
1998) of the form 

SFR = ( 5 mf x 10“ 10 L ir (1) 

where SFR is in units of M 0 yf' and Lir is in units of Lq. 
The normalization <5 mf is sensitive to assumptions about the 
initial mass function, the duration of star formation, and dust, 
discussed in more detail in §4.2. However, it is generally cal¬ 
culated to be a factor of order unity (0.8 < dvir < 2 in, e.g. 
Scoville & Young 1983; Thronson & Telesco 1986; Kenni¬ 
cutt 1998; Barger et al. 2000; Rowan-Robinson 2000; Omont 
et al. 2001). 

We adopt a Chabrier (2003) initial mass function, which 
conventionally entails <5 mf =1.0 (e.g. Magnelli et al. 2012; 
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Carilli & Walter 2013), as was used by Behroozi et al. 
(2013a). 

To our knowledge, there has not been a recent compre¬ 
hensive study of the expected range of <5mf for high-redshift 
galaxies, and such a study is beyond the scope of this pa¬ 
per. In order to remain consistent with the values quoted 
above, we adopt a log-normal prior of log^F = 0.0 ± 0.3 
(5 mf ~ 1.0^-g). Note that this means the prior’s ±3cr interval 
spans nearly 2 dex. 

Eq. 1 could also be written as SFR = <5mf^sb x 10 _10 Lir, 
where 5 sb explicitly accounts for the fraction of Lir due 
to starbursts, as opposed to active galactic nuclei or other 
sources. For simplicity, we have absorbed <5 sb into Smf* and 
any scatter about Eq. I due to variation in galaxy type is as¬ 
sumed to be absorbed into ctsfr (§2.2.2). 


to Lqo (units of L 0 ) is 


Leo = 4.9 x 10 5 Lq 


^CO. rest 

115.27 GHz 


3 


( ^ 


'CO 


V Kkm s 


-l 



(4) 


where ^co. rest = 115.27 GHz is the rest-frame frequency of 
the CO transition. 

To resummarize the model: 


1. Halos —> SFR: Get SFR(M,z) from the results of 
Behroozi et al. (2013a) 

2. Add log-scatter, ctsfr 

3. SFR ->• L ir : Get Lir from SFR = <5 M f x 10- 10 L ir 

4. Lir L' co : Get L' co from log Lir = a log L' co + /3 


2.2.4. Infrared Luminosity and CO Luminosity 

To convert infrared luminosity to CO luminosity, we as¬ 
sume a power-law relation of the form 

log Lir = a log L' co + ft (2) 

where L JR is in units of L 0 , and L' co is in units of Kkm s -1 pc 2 
(observer units for velocity- and area-integrated brightness 
temperature). 

For our fiducial model, we use the fit from Carilli & Walter 
(2013), which found from a census of high-redshift galaxies 
a = 1,37(±0.04) and /? = -1,74(±0.40). 

More generally, however, this relation is a rough proxy 
for the relation between star formation and molecular gas, 
and has been fit to data from high-redshift galaxies (z > 1, 
but typically few if any galaxies beyond z ~ 3) in a num¬ 
ber of studies. These studies have found, e.g., (a,/if) val¬ 
ues of (1.13,0.53), (1.37,-1.74), (1.00,2.00), (1.17,0.28) 
(Daddi et al. 2010; Carilli & Walter 2013; Greve et al. 2014; 
Dessauges-Zavadsky et al. 2015, respectively). These values 
are closely fit by the line 

a ~ 0.10/3+1.19 (3) 

as they are generally derived from similar samples of the 
most luminous, and therefore detectable, high-- galaxies. We 
choose deliberately loose priors of a = 1.17 ±0.37 and j3 = 
0.21 ±3.74 by taking the mean of those four values and using 
their full range as the ler spread. 5 

We also add log-normal scatter to Leo, fiducially parame¬ 
terized by <7i co = 0.3. This is effectively the scatter in Leo for 
a given L JR or SFR (they scale linearly and so are separated in 
log-space by a constant). However, equal values of CT£ co and 
ctsfr do not necessarily have the same effect on the signal, 
since SFR and Leo are not required to scale linearly. In prin¬ 
ciple, though, f t/^q and ctsfr could be combined into a single 
scatter parameter. We choose not to do so because empirical 
constraints, which inform our priors, are generally on the two 
separate parameters. 

The prior on this parameter is ctl co = 0.3 ±0.1, with the 
requirement that ctl co > 0. This is consistent with the results 
of aforementioned studies, where a scatter has been quoted. 

Note that the conversion from L' co (units of K km s -1 pc 2 ) 

5 While this paper was under review, the values of (a,/3) in Dessauges- 
Zavadsky et al. (2015) were updated from (1.19,0.05) to (1.17,0.28). We 
have updated text and figures in this work where appropriate. 


5. Add log-scatter, oq co 

with fiducial parameter values: 

(’’sfr = 0.3, = 0.3, 

S M f= 1.0, a = 1.37, /? = -1.74 

Figure 2 shows the combined result of these steps, plotting 
the mean Lco(4L/,) relation from our fiducial model, as well 
as the equivalent relation from previous studies. Notably, Leo 
in this model is not linear in M, a simplifying assumption that 
has been adopted in some previous studies. This is a con¬ 
sequence of our choice to model the average population in¬ 
cluding quiescent galaxies, which results in a flatter function 
of mass at the high-mass end than a model constrained us¬ 
ing only detected star-forming galaxies. This can affect the 
shape of the power spectrum, even if the mean CO brightness 
temperature is the same. 

Our model may be seen as combining the two models from 
Pullen et al. (2013), “A” and “5”, in that we still start with 
dark matter halos (their model A) but calculate empirically 
constrained SFRs (their model B ). In short, we combine the 
machinery for A with the motivation for B. 

2.2.5. Generating Intensity Maps and Power Spectra 

Once halo CO luminosities are calculated, we generate 3D 
intensity maps (spectral data cubes) and power spectra. We 
calculate intensities as brightness temperatures, since we ex¬ 
pect observation and calibration to be in terms of that quantity. 

To do this, the halos are binned by their positions on an 
RA x Dec x r/obs grid, where r^bs = r'co, rest /(1 +z). From the 
total CO luminosity in each cell, brightness temperatures are 
calculated using the Rayleigh-Jeans relation 


Here, /, / ohs = Lco/4-KDjd,^ where D L is the luminosity dis¬ 
tance to the source, and is the width of the frequency 
channel (spectral resolution element), specified in §2.4. This 
assumes that the CO line profile is approximately a delta 
function. In this study, we have not included the effects 
of Doppler broadening and redshift-space distortions, which 
would somewhat alter the line-of-sight signal, but they are left 
for future studies. 

We calculate power spectra P( k) by converting the grid to 
3D comoving units and squaring the discrete Fourier trans¬ 
form of the temperature cube. More details are provided in 
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Appendix B. Averaging the 3D power spectrum P(k) in radial 
bins gives us the spherically averaged power spectrum, P(k). 
In this paper, we plot on the quantity A 2 (k) = k s P(k) flrf in¬ 
stead of P(k) directly, mainly to compare with previous stud¬ 
ies. The physical meaning of A 2 (k) is the amount of variance 
in T contributed per logarithmic interval in k. 6 

We assume white (Gaussian) instrumental noise in these 
simulations. As with all measurements, actual data will have 
additional non-ideal signal components, but we do not at¬ 
tempt to add these instrument-specific terms here. By us¬ 
ing P(k) as calculated above, we implicitly assume that the 
CO power spectrum has been cleaned of instrument noise. In 
this case, we can imagine having subtracted a constant noise 
power spectrum P n (k) oc a 7 (Eq. C2) from the noisy spectrum, 
leaving only residual fluctuations from noise variance. Error 
bars on P(k) have been calculated to account for this noise 
variance, as well as resolution limits and sample variance of 
P(k). We detail this calculation in Appendix C. 

When showing predictions for the power spectrum, we av¬ 
erage over 100 realizations of the same volume, where each 
realization is taken from a separately generated lightcone. 
This averaging smoothes out fluctuations in the power spec¬ 
trum that are purely sample variance. While those fluctua¬ 
tions due to sample variance are important to consider, they 
distract from the average power spectrum predicted by the 
model, which is what we aim to show unless otherwise speci¬ 
fied. 

2.3. Part 2: MCMC Inference From Mock Signal 

Modeling the expected signal is a useful first step. How¬ 
ever, given the significant empirical uncertainties of this (or 
any) model, there is a limit to how useful a single prediction 
can be. 

A question of at least equal importance is how intensity 
mapping can inform our understanding of the underlying 
galaxy population. Additionally, it is worth addressing the 
implications for galaxy properties if no clear CO signal is de¬ 
tected by the experiment: given the uncertainties, this is a pos¬ 
sibility. While a non-detection might not be the most desired 
outcome, we wish to know whether it would be informative 
and therefore scientifically interesting. To that end, we pro¬ 
duce the following three “mock” power spectra: 

1. P(k) from our fiducial L C q(M) model. 

2. P(k )ss 0 within random fluctuations ( o P ~ P n /^/N modes; the 
second term of Eq. Cl). This represents a signal consistent 
with a non-detection after subtracting a flat noise spectrum. 

3. P(k) from the simulated data of Obreschkow et al. (2009b), 
which were generated with a model different from ours. 

Combined with priors informed by existing observations, we 
use the intensity mapping power spectrum to infer constraints 
on our model parameter space for all three of these models. 

We generate these three mock spectra and their error bars 
for both proposed instruments discussed in the following sec¬ 
tion (§2.4), for a total of six hypothetical signals. The purpose 
of the fiducial signal is to see how closely we can recover the 
“true” parameter values. The purpose of the “zero” power 

6 A-(tr) is sometimes called the “dimensionless” power spectrum if the 
real-space quantity being analyzed is unitless, e.g. number overdensity. We 
avoid this label here because the real-space quantity of interest is brightness 
temperature, so A 2 (k) has units of /.iK-. 


spectrum is to investigate what constraints a non-detection 
places on the galaxy population. The purpose of the signal 
from Obreschkow et al. (2009b) data is to check whether our 
results are still reasonable when an independent model is used 
to generate the CO signal. 

We use emcee 7 , an implementation of an affine-invariant 
ensemble Markov chain Monte Carlo (MCMC) algorithm 
(Foreman-Mackey et al. 2013; Goodman & Weare 2010), 
to sample the posterior distributions in our model parameter 
space. We provide more details, including the likelihood, in 
Appendix B.3. 

A summary of the model parameters and their priors is pro¬ 
vided in Table 1 . All parameters are sampled in linear space 
except for <5 mf, which is sampled in log space. 

In this study we model osfr and cr Lco as uncorrelated, so 
from the steps described in §2.2, one could write the total 
scatter as er tot = (osfr / 0(2 + a L co ^^ 2 ■ However, if we consider 
perfectly correlated or anticorrelated scatter, the total scatter 
would instead be a tot = (tsm/oi+a^ or cr to t = IcrsFR/a-crzcol. 
This would change the fiducial value (cr tot ~ 0.37) to cr tot ~ 
0.52 or crtot ~ 0.08, respectively. Insofar as we only care about 
the total scatter, these extreme cases should be adequately 
covered by our priors. Strongly correlated scatter could plau¬ 
sibly affect the inferred values of other parameters, but we do 
not model this for the present study. 

In this work, we have chosen to err on the side of broad pri¬ 
ors. A more thorough analysis of the uncertainties, system- 
atics, and details outside the scope of this study is warranted 
and welcome in future papers. 

2.4. Instrument and Survey Design 

To connect this work to proposed and feasible observations, 
we consider two examples of possible dedicated instruments 
for CO intensity mapping that can be built with current detec¬ 
tor technology. In our calculations, the instrument and survey 
parameters determine the size and resolution of our 3D inten¬ 
sity map, as well as the error bars on the power spectrum. 

An essential requirement of intensity mapping is high fre¬ 
quency resolution. In contrast to current measurements of ex- 
tragalactic backgrounds—such as the cosmic infrared back¬ 
ground (e.g. with Spitzer or Herschel ) or the cosmic mi¬ 
crowave background (CMB)—which characterize excess ra¬ 
diation in relatively broad bands, the redshifted CO inten¬ 
sity from galaxies will have narrow fluctuations in frequency 
space, reflecting cosmological structure along the line of 
sight. For reference, for CO(l-O) emission from z, = 2.5 
(r'obs ~ 33 GHz), a spatial depth of 1 comoving Mpc corre¬ 
sponds to a frequency width of 8.25 MHz, which requires 
high-resolution spectroscopy to resolve. 

Beyond that, both instruments we consider are single-dish, 
ground-based telescopes covering a band of 30-34 GHz in the 
Ka microwave band. These frequencies correspond to red- 
shifted CO(l-O) emission from z « 2.4-2.8. 

The first instrument is envisioned as a pathfinder exper¬ 
iment, with parameters similar to the CO Mapping Array 
Pathfinder (COMAP), currently under development. Its main 
goal would be to detect or constrain the CO signal from galax¬ 
ies at z ~ 2-3. We have assumed an aperture of ~6 m with 
a frequency resolution of R = v 0 \, s /8 v ~ 800, with a 2.5 deg 2 
field of view. For the rest of the paper, we refer to this instru¬ 
ment as the “pathfinder” experiment. 

7 http : //dan . iel . fm/emcee 
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an intensity mapping experiment. There are at least a few dis¬ 
tinct ways to address this: 

• Mean CO intensity vs instrument noise. One can com¬ 
pare the mean CO brightness temperature {T co ) to the in¬ 
strument noise fluctuations cr n (Eq. C6), which is a useful 
first check (e.g. Carilli 2011, at z > 6). However, this is not 
exactly the signal being measured, and it is possible to have 
an instrument where (Tqo) < er„ in individual map “vox¬ 
els” (as with the pathfinder experiment in this work), and 
still measure cosmological intensity fluctuations on larger 
scales. On those scales (across multiple voxels), random 
noise fluctuations should average out, while cosmological 
fluctuations should remain. Moreover, an actual observa¬ 
tion would require subtraction of continuum foregrounds, 
making a direct measurement of (Tqo) (the k = 0 mode) dif¬ 
ficult. 

• Measuring spatial fluctuations in the CO signal. In fact, 
we are interested in detecting the spatial structure of the 
signal. Resolved 3D spatial structure is essential for inten¬ 
sity mapping. Without it, there is no way of determining 
whether a signal is CO from high-redshift galaxies, which 
trace cosmological structure along the line of sight, or sim¬ 
ply an unwanted systematic (e.g. synchrotron foregrounds, 
which are generally smooth in frequency space). For an ex¬ 
periment with the goal of making a detection, we can sum 
in quadrature the detection significances of each measured 
mode (Pullen et al. 2013; Breysse et al. 2014, Eq. 7). 


Fig. 3.— Fiducial CO power spectrum and detection significance. Top: 
Spherically averaged CO power spectrum from our fiducial model (§2.2), as 
measured by the pathfinder (red) and full (blue) experiments and averaged 
over four separate sky patches. Solid lines show the measured average power 
spectrum. Shaded regions show the lcr uncertainty from sample variance. 
Dashed lines indicate the lcr limits from thermal instrument noise. The up¬ 
turn in these limits at high k (small scales) arises from resolution limits (finite 
beam and channel width). Bottom: Detection significance of the power spec¬ 
trum as a function of k. The total detection significance (Eq. 7) is 7.89 (144) 
for the pathfinder (full) experiment. 

The second instrument is a more advanced experiment, 
which would aim to measure the CO signal with greater sen¬ 
sitivity, over a wider area, and with better resolution. To this 
end, we have assumed an aperture of ~12 m aperture with 
R ~ 3200, with a 6.25 deg 2 field of view. For the rest of the 
paper, we refer to this instrument as the “full” experiment. 

The relevant parameters from the two instruments are sum¬ 
marized in Table 2. For the purposes of our calculations, we 
assume that each instrument will pursue a survey that ob¬ 
serves four separate patches on the sky, a necessary survey 
strategy due to the limitations of ground-based observations. 
However, we can combine the power spectra from each patch, 
reducing the error bars on the power spectrum (see Appendix 
C) by a factor of approximately \/4 = 2. 

Throughout this paper, we use units of brightness tempera¬ 
ture. If desired, the unit conversion from brightness tempera¬ 
ture T to flux density S can be written as 


S = 0.108 mJy beam 1 


^obs 

32 GHz 


2 


^FWHM 



2.4.1. Sensitivity 

We have given specific parameters for a CO intensity map¬ 
ping experiment, but in this section we discuss, more broadly, 
the question of how we should characterize the sensitivity of 


• Usefulness for inferred constraints. This requires a prop¬ 
erly targeted range of k. To measure cosmological large- 
scale structure, e.g. baryon acoustic oscillations or the lin¬ 
ear power spectrum, one would like to cover large areas, 
focusing mainly on the clustering modes at low k. To make 
inferences about the galaxy population, one also needs sen¬ 
sitivity to the shot noise modes at high k. However, the ex¬ 
act transition between the clustering and shot noise scales 
depends on the details of the underlying galaxy population, 
as shown in §3.2. One needs sensitivity to P(k) in both 
regimes, either from a single experiment or from multiple 
experiments targeting different scales. Information about 
the underlying galaxy population is contained in the rela¬ 
tive strength of the power spectrum on different scales. 

3. RESULTS 

3.1. Fiducial Prediction and Detection Significance 

Figure 3 shows the fiducial power spectra obtained from 
our mock intensity maps (e.g. Fig. 1), as well as the mode 
detection significance as a function of k. As has been noted 
by previous studies, the power spectrum may be thought of, in 
the context of the “halo model” formalism (Cooray & Sheth 
2002), as having a clustering component (at low k ) and a shot 
noise component (at high k). 

In averaging P(k), we choose bin widths of A k = 2ir/L m i n , 
where L rnln is the shortest comoving dimension of the survey 
volume, specifically A k = 0.038 Mpc -1 (for the pathfinder) 
and 0.024 Mpc -1 (for the full experiment). Note that the 
apparent deficit in the pathfinder’s low-k power is because 
its wider, linearly spaced k-space bins cause the the average 
power in the lowest bin to appear underestimated. 

Averaging the power spectrum over four identical patches, 
as mentioned in 2.4, we find that the maximum detection sig¬ 
nificance P(k)/crp(k) for any single k mode is 4.46 and 35.2 
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Fig. 4.— Range of power spectra spanned by ±1 cr uncertainties on the 
mean SFR (M,z) from Behroozi et al. (2013a). At a fixed scale k, these span 
power spectra values of approximately d=0.2 dex. No halo-to-halo scatter in 
Leo has been included here (the effect of scatter is shown in Fig. 5). 

for the pathfinder and full experiments, respectively. How¬ 
ever, we can quantify the total detection significance of the 
power spectrum as (e.g. Pullen et al. 2013; Breysse et al. 
2014) 

snrL = E 

i 

where k, are the discrete values of k at which the power spec¬ 
trum is calculated, separated by A k. 

Assuming uncorrelated errors, we find the total signal-to- 
noise ratio (SNR) to be 7.89 and 144 for the pathfinder and 
full experiments, respectively. However, those values are only 
for our specific choice of modeling and observing parameters. 
In general, SNR tot may be improved by increasing bandwidth 
or by decreasing noise fluctuations (the latter may be achieved 
by increasing observing time, increasing the number of feeds, 
or decreasing system temperature). 

It is also possible to optimize the survey area to maximize 
SNRtot (see, e.g. Breysse et al. 2014), though note that the 
exact optimal area depends on P(k ) and is therefore model- 
dependent. When we perform this calculation, we obtain 0.6 
and 9.2 deg 2 for the pathfinder and full experiments, respec¬ 
tively. These values differ from our choices of of 2.5 and 6.25 
deg 2 , but adopting the optimal values does not yield dramatic 
improvements (increasing SNR tot to 9.48 and 146). Through¬ 
out this study, then, we maintain the pathfinder and full exper¬ 
iment parameters in Table 2. See Appendix D for more details 
and considerations. 

3.2. Varying Model Parameters 

The uncertainties in SFR alone allow the power spectrum 
to span ±0.2 dex in amplitude. This can be seen in Figure 
4, which shows the range of power spectra expected from the 
the lu posterior on SFR(M/,,z) (Behroozi et al. 2013a). As 


noted, the width of this interval is fairly consistent in log space 
over all halo masses at the redshifts of interest (~ ±0.15 dex). 
As a result, varying SFR(7T7/,.z) within the posterior can be 
approximated as a simple rescaling, and the effect of doing 
so is degenerate with that of 5 mf- In the MCMC inference 
procedure (§2.3), we have absorbed the effect of this rescaling 
entirely into 5 MF . 

Figure 5 shows the effect of varying halo-to-halo scatter via 
(Jsfr and oy, co . Increasing scatter increases the relative am¬ 
plitude of the shot noise component of the power spectrum, 
compared to the clustering “bump” at low k. In the limit of 
high scatter, P(k) ~ const, or A 2 {k) oc k 3 . Notably, halo-to- 
halo scatter introduces sample variance in the signal within 
the same volume. However, this effect only appears to be sig¬ 
nificant for very high scatter: cts FR ~ 1.0 dex or ay, co ~ 1.0. 
These values are not expected to be realistic since they imply 
a 68% scatter of 2 dex. 

We do not show the effect of varying the SFR-IR normal¬ 
ization <5 mf, since it is simply a scaling factor: with other pa¬ 
rameters fixed, a higher value of <5 M f results in fainter CO 
brightness. For a given galaxy, a high value for (5 MF (kP 1) 
would mean a large amount of star formation per IR luminos¬ 
ity. This could could indicate very low dust content or very 
short star formation timescales. 

Figures 6 shows the effect of varying a and (3 in the ±i R - 
L' c0 relation. Rewriting Eq. 2 as logL) :0 = o r 'tlog L| R — /j) 
makes the effect of each parameter more apparent. Indepen¬ 
dently increasing a or /? decreases L' co for a given L IR . How¬ 
ever, increasing a also weights halos with low ±i R (~ low 
SFR) more relative to halos with high ±i R (~ high SFR). De¬ 
pending on which halos are occupied by high-SFR galaxies, 
this can affect the shape of the power spectrum, particularly 
the clustering signal. This is a notable difference that can arise 
in predictions of power spectra if SFR is not modeled as sim¬ 
ply linear in M/, (in our model, the average SFR turns down¬ 
ward just above M/, ~ 10 12 M 0 , a consequence of the now 
well-established downturn in the ratio of stellar mass to halo 
mass at a given halo mass). 

3.3. MCMC Inference 

Figure 7 shows the posterior distributions for all model pa¬ 
rameters, as inferred from the “fiducial” and “non-detection” 
power spectra. In general, the full experiment places stronger 
constraints on the model parameters than the pathfinder. 

Figure 8 shows constraints on the Li R -L' co relation and the 
CO luminosity function, from both the pathfinder and full ex¬ 
periment observing our fiducial signal. 

In the Li R -Lq 0 plots, we have overplotted the four empiri¬ 
cal fits mentioned in §2.2.4 (Daddi et al. 2010; Carilli & Wal¬ 
ter 2013; Greve et al. 2014; Dessauges-Zavadsky et al. 2015). 
The inferred relation is consistent with the fiducial line (Car¬ 
illi & Walter 2013) within ler, and consistent with all four 
lines within 2a. The notable point here is that while the four 
empirical fits were extrapolated from a limited, bright sam¬ 
ple of galaxies (hence their convergence around ~ 10 9 - 10 u) 
K km s _1 pc 2 ), the relation inferred from the intensity map 
was constrained by probing faint populations directly , albeit 
in integrated emission. 

In the plots of CO luminosity function (LF), we recover 
the “true” CO luminosity function (underlying the fiducial 
intensity map) to within ler over nearly all luminosities. 
Here, the full experiment produces tighter constraints than the 
pathfinder experiment, as expected. Note that the “true” LF 
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Fig. 5.— Effect on the CO power spectrum of halo-to-halo scatter as parameterized by: SFR given A//, (ctsfr> left) or Eco given SFR (c tl cq , right). At 
csfr^Lco « 1 dex and higher, two things are evident: (1) the power spectrum begins to look like a pure shot noise spectrum, since any clustering signature is 
increasingly buried by the large halo-to-halo scatter, and (2) the scatter introduces significant variance into the power spectrum. This variance is demonstrated 
for the a = 1.0 case by plotting 100 individual power spectra in the top plot, as well as their shaded 95% interval in the bottom plot. However, this scenario is 
probably extreme and unrealistic, because it implies a =blcr scatter over two orders of magnitude. In this plot, the labeled lines are the mean of 100 power spectra 
from the exact same halos (rather than multiple realizations of the survey volume) to isolate the variance introduced by halo-to-halo scatter alone. 



Fig. 6.— Power spectra due to varying Lir-Lco relations. Solid colored 
lines show literature values of a, (3 (Daddi et al. 2010; Carilli & Walter 2013; 
Greve et al. 2014; Dessauges-Zavadsky et al. 2015). Gray dashed lines show 
the power spectra from a = 0.5,1.0,1.5,2.0 (dashed lines) and corresponding 
P assuming a and /3 covary as in Eq. 3. As a increases, smaller and fainter 
(in IR) galaxies contribute a greater proportion of the total CO luminosity, 
and the relative strength of the low-k clustering signal increases. 

lies outside of the ler interval at the bright end. While this is 
not a very strong deviation, we expect that intensity mapping 
will not be especially sensitive to the bright end of the LF if 


somewhat fainter populations dominate the overall signal, as 
they do here. 

We also show 2er posterior constraints on the CO luminos¬ 
ity function in the event that no clear signal is detected (Fig. 8, 
bottom right). While these constraints visually resemble up¬ 
per limits, they are actually allowed regions and are therefore 
sensitive to our choice of priors. We include further discus¬ 
sion of these priors in §4.3. 

To compare the FF constraints with a hypothetical blind 
survey for CO in individual galaxies, we have marked the 
faintest CO luminosity that could be detected by the VLA 
in 1 and 100 hr of on-source time. For the calculation, we 
used the VLA exposure calculator 8 , assuming an observing 
frequency of 32 GHz (z « 2.6), 100 km s -1 channel, and 5er 
detection threshold. 100 hours of on-source time can probe 
below the “knee” of the luminosity function predicted by this 
model. However, the relatively small field of view (~ 1.3' pri¬ 
mary beam at 32 GHz) means that sample variance in such a 
blind survey will be significant. 

3.3.1. Cross-check with Mock Data from Obreschkow et al. ( 2009) 

The model used to generate the mock signal is the same 
model used in the MCMC procedure. It is reassuring, but per¬ 
haps unsurprising, that we recover the “true” CO luminosity 
function when fitting this model. What if we use a population 
of CO galaxies that was generated from an independent model 
applied to a different simulation? 

Various studies have provided physically motivated model¬ 
ing of Leo at high redshift, whether through direct hydrody¬ 
namic simulation or semi-empirical modeling (e.g. Popping 
et al. 2014; Lagos et al. 2015). For this exercise, meant to pro- 

8 https://obs.vla.nrao . edu/ect 
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Fig. 7.— MCMC posterior distributions on the parameter space of our model. Black contours show prior distributions. Red contours show constraints from the 
pathfinder experiment, while blue contours show constraints from the full experiment. Crosshairs indicate the values of our fiducial model. Marked points in the 
a-(3 plot (bottom row, second from right) show values obtained by previous studies (see §2.2.4). Left: Posterior distributions inferred from our fiducial signal, 
as observed by both the pathfinder (red) and full (blue) experiments. Right: Posterior distributions inferred from a non-detection by the same two experiments. 


vide a sanity check rather than a prediction, we assume the un¬ 
derlying galaxy Leo’s are given by the results of Obreschkow 
et al. (2009a,b), who use a semi-analytic model to predict 
galaxy CO luminosities from dark matter halos. From their 
publicly downloadable lightcone data 9 , we generate CO inten¬ 
sity maps and rerun the same analysis, for both the pathfinder 
and full experiments. See Appendix E for more details. 

Figure 9 compares the actual CO luminosity function in the 
volume with the inferred luminosity function from intensity 
mapping. The weaker signal is consistent with non-detection 
for the pathfinder. However, it is detected by the full experi¬ 
ment, and within the 95% credible intervals, the inferred LF 
is consistent with the "true" LF, while our original model’s 
LF is largely excluded. While this is not an exhaustive check, 
it is an indication that our model appears flexible enough to 
accommodate a reasonable range of predicted LFs. 

4. DISCUSSION 

4.1. Relation to Previous Work 

Our work here naturally extends the various studies men¬ 
tioned in §2. 1 . In modeling, our main improvements are (1) an 
empirically constrained connection between halos and mean 
SFRs and (2) an Lir-Lco connection that is fit to a larger set 
of observations. These should model more accurately the con¬ 
nection between underlying large-scale structure and the pre¬ 
dicted CO intensity map. Qualitatively, the shape of the over¬ 
all Lco(M/,) relation is somewhat similar to that of Gong et al. 
( 2011 ), though their focus was on z > 6 . 

Figure 10 shows the power spectrum of our fiducial model 
compared to other models at similar redshifts (see Fig. 2). 
Note that the power spectra of those models were calculated 
with a halo mass cutoff of Mco, min = 10 9 M Q , lower than our 
fiducial 1 O 1 () M 0 , to better compare with their original pre- 

9 http : //s-cubed . physics . ox . ac . uk/s3_sax 


dictions. However, this does not affect the main qualitative 
point of the plot: the nonlinear form of our Lco(M /,) affects 
the shape of the power spectrum — in this case, increasing 
the relative power in clustering (low-k) modes. 

The novel aspect of this work is an attempt to use CO inten¬ 
sity mapping — or rather, simulated observations thereof — 
to make inferences about the underlying galaxy population. 
This is particularly significant because, despite the impres¬ 
sive capabilities of ALMA and VLA, it will still be difficult 
for them to probe complete samples of the faint end of the CO 
luminosity function. Considering that CO is the most accessi¬ 
ble tracer of molecular gas, relying solely on galaxy surveys 
may limit our ability to trace the gas-SFR connection at high 
redshifts. 

In our approach, we have used full /V-body simulation, 
which is more computationally expensive than analytically 
approximating the power spectrum, even if existing halo cat¬ 
alogs for a simulation are already available. Indeed, with a 
simple Leo oc M model, the power spectrum can be quickly 
and analytically calculated from the halo model formalism 
(Cooray & Sheth 2002). This is how the Leo oc M power 
spectra (i.e. previous models) of Figure 10 were calculated. 
However, directly simulating these halos is a method that 
more naturally accommodates information about galaxy en¬ 
vironment and merger histories. For example, merger-driven 
starbursts could be identified and modeled separately from the 
“normal” star-forming population. Modeling of that nature is 
left for future studies. Galaxy formation involves many com¬ 
plex processes with diverse observational signatures, and this 
approach provides a framework for all of those facets to be 
studied more directly. 

4.2. Modeling Uncertainties 

One goal of intensity mapping is to probe high-redshift 
galaxies too faint to be individually observed in large num- 
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Fig. 8.— The properties of galaxy populations, as inferred by the pathfinder (red) and full (blue) experiments. Top left: Inferred I a constraints on the mean 
L]u C;o relation. Note that while the pathfinder may appear to yield comparable constraints to the full experiment, this is an artifact of our somewhat simple 
priors and likelihood (see §3.3). The parameters a and ft are individually more constrained by the full experiment (see Fig 7). For comparison, we also plot 
derived fits from Daddi et al. (2010), Greve et al. (2014), and Dessauges-Zavadsky et al. (2015). Top right: 95% posterior constraints for the Lir- L' cq relation 
from a non-detection. Bottom left: Inferred CO luminosity functions. Colored intervals show the median and 95% credible interval for inferred CO luminosity 
functions. For comparison, we show observational constraints at z ~ 2.75 from a blind scan of the Hubble Deep Field (Walter et al. 2014), as well as the model 
predictions of Obreschkow et al. (2009a) at similar redshifts. The vertical dashed lines marks the VLA 5cr detection threshold assuming 1 and 100 hr of on-source 
time and a 100 km s -1 channel width at 32 GHz. Bottom right: 95% posterior constraints on the CO luminosity function from a non-detection. 


bers. By nature, then, these are currently poorly studied pop¬ 
ulations, and any attempt to empirically model their Leo’s 
must rely on data from (1) galaxies thought to be low-redshift 
analogs, e.g. local dwarfs and spirals; (2) observed galaxies 
at high redshift, which may be atypically luminous or star¬ 
forming; or (3) galaxies that might be “typical” at high red- 
shift, but in small and incomplete samples. 

Predictions about such galaxies are subject to significant 
uncertainties about the formation and evolution of high- 
redshift galaxies. In particular, the nature of the ISM at high 
redshift is not fully understood. In the very early universe, one 
might expect significantly lower metallicity and dust shielding 
in the ISM, resulting in lower CO abundance. However, in ac¬ 
tive star-forming regions, young massive stars should quickly 
pollute their surroundings with metals. These various effects 
have not been well quantified in typical high-redshift galax¬ 
ies and are just beginning to be probed by individual galaxy 
observations with, e.g., ALMA and VLA. 


In our model, L IR correlates with SFR (Eq. 1). This is be¬ 
cause we need an empirical tracer of star formation to link 
halo SFR with CO luminosity. Using Lir assumes that star 
formation occurs in the presence of dust, and that Lir is ther¬ 
mally radiated from dust heated by massive, young stars. 

The normalization <5 mf, as well, is sensitive to star forma¬ 
tion history and the initial mass function, both of which may 
be different at high redshift. It should also be noted that 
a tracer like Lir may overestimate the instantaneous SFR if 
star formation occurs in short “bursty” episodes (<C 100 Myr, 
e.g. Dominguez et al. 2015) or in recently quenched galaxies 
(Hayward et al. 2014). 

Additionally, although we have taken L !R to be the total lu¬ 
minosity in the 8 —1000 /im band, note that this is not always 
the directly observed quantity. Lir may be converted from a 
measurement of flux in a narrower IR band. Any systematic 
errors introduced by this fact will probably not be dominant, 
but it is something to remember going forward. 


4.2.1. Using Lir to trace SFR and Leo 


4.2.2. CO Line Luminosities 
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Fig. 9.— CO luminosity function from the data of Obreschkow et al. 
(2009b), as inferred by the pathfinder and full experiments. Our fiducial 
model has been plotted for comparison, but note that the inferred LF matches 
the Obreschkow et al. (2009b) curve more closely now while excluding our 
fiducial model (as should be expected). For more details, see Fig. 8. 



Fig. 10.— Comparison of the fiducial power spectrum in this work with that 
of previous models (Righi et al. 2008; Visbal & Loeb 2010; Pullen et al. 2013, 
see Fig. 2), which either modeled Lco c* or have been approximated as 
such here. Our model does not have a linear Lco(Mh), and this qualitatively 
affects the shape of the resulting power spectrum. Nevertheless, our fiducial 
prediction appears to be consistent with the range of previous predictions. 

Existing high-redshift data for CO(l-O) luminosity are also 
not free of uncertainty. Some values for the CO(l-O) lumi¬ 
nosity were not directly measured, but instead inferred from 
higher-order transitions that were accessible in an available 
observing band (e.g. Genzel et al. 2010; Tacconi et al. 2013). 
Inferring CO(l-O) luminosities from such observations re¬ 
lies on assumptions about spectral line energy distribution 
(SLED), i.e. the relative luminosities of each CO line. In 
the optically thick, high-temperature limit, all lines have the 
same brightness temperature, or equivalently luminosities of 
Lco(j^j- 1) oc J 2 Lcoii-o), but this generally overestimates 


high-/ luminosities relative to Lco(t-O)- 

While some recent work has focused on characterizing the 
SLED across a sample of galaxies (e.g. Greve et al. 2014; 
Narayanan & Krumholz 2014), the scatter is large and predict¬ 
ing the SLED in an individual galaxy is not straightforward. 
On the other hand, it may be possible to constrain the average 
SLED by cross-correlating two or more intensity maps that 
target the same cosmological volume through different CO 
transitions. 

4.3. Inferred Constraints on Galaxy Populations 

The Lir-L' co relation (Eq. 2) has been fit to observations 
in previous literature (Daddi et al. 2010; Dessauges-Zavadsky 
et al. 2015; Carilli & Walter 2013; Greve et al. 2014). In Fig¬ 
ure 8, we have plotted this relation, both as derived in the 
literature and as inferred from a hypothetical CO signal. To 
the extent that the L tR -L' C0 relation traces the star formation- 
molecular gas relation, it is possible for intensity mapping to 
characterize the latter in unobserved galaxy populations. 

Our results also constrain the CO luminosity function. Cur¬ 
rently, this is not well characterized at the redshifts in this 
study because of the long integration times required to detect 
CO in a single galaxy. We have overplotted the results of Wal¬ 
ter et al. (2014), which estimated the CO luminosity function 
from a blind scan for CO lines in the Hubble Deep Field. 

The noteworthy point here is that these constraints, es¬ 
pecially for the faintest galaxies, arise from directly imag¬ 
ing their emission, albeit in aggregate. By contrast, because 
galaxy surveys preferentially probe the brightest galaxies, in¬ 
ferences about the faintest galaxies need to be made by, e.g., 
extrapolation. With intensity mapping, part of the informa¬ 
tion we receive comes directly from the (integrated) faint end 
of the luminosity function. 

In future analyses, one can certainly consider tighter priors 
than we have chosen here. In the event of a detection of an 
intensity mapping signal, we expect narrower priors to yield 
narrower constraints on the CO luminosity function. In the 
event of a non-detection, we would actually expect narrower 
priors to yield less stringent constraints, i.e. a higher upper 
bound on the allowed luminosity function. This is because 
our broad priors on a and j3 (of the L ir -L C q relation) actu¬ 
ally allow a significant region of parameter space (high a and 
3) that would be in significant tension with existing Lir-Lco 
data, even though it yields undetectable CO power spectra, 
consistent with the “observed” intensity map. 10 As a result, 
given a mock non-detection, a region of parameter space pro¬ 
ducing unrealistically faint populations is still allowed in the 
posterior region. Tighter priors would exclude that unrealistic 
parameter space, raising the upper bound on the allowed lu¬ 
minosity function. All the same, the plots in Fig. 8 do demon¬ 
strate that a non-detection of an intensity mapping signal can 
constrain the faint end of the luminosity function. Alterna¬ 
tively, a more complete likelihood calculation—for example, 
simultaneously fitting the L ]K -L' CQ relation to galaxy survey 
and intensity mapping data—should also yield stronger and 
more robust constraints. 

One parameter we have not allowed to vary is Mqo, m i n , 
the minimum mass of CO-luminous halos, which we fix at 
10 10 Mg. In principle, this could also affect the signal. How¬ 
ever, as we note in Appendix A, lowering the mass cut- 

10 This also appeal's to be why the Ljr-LC constraints in Fig. 8 do not 
seem to be improved by the full experiment, a somewhat counterintuitive 
result. 
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off yields diminishing returns in our fiducial model because 
Lco(M) decreases toward low M faster than the halo number 
density increases. Any variation in the signal we might get 
from decreasing Mqo. m in below below 10 10 M 0 is expected 
to be eclipsed by other uncertainties in the model. However, 
one can certainly imagine more realistic models than those 
using an abrupt halo mass cutoff. We leave investigation of 
these models to future studies. 

A consideration in making meaningful inferences is that the 
spherically averaged power spectrum essentially carries two 
independent sources of information: (1) the amplitude of the 
clustering (low-A') signal, and (2) the amplitude of the shot 
noise (high-k) signal. This restricts the complexity of mod¬ 
els we can expect to meaningfully constrain. However, there 
are other effects that are not modeled here but will be imple¬ 
mented as refinements in a future study. These include red- 
shift space distortions and line broadening, which will distort 
the signal along the line of sight. This slightly alters the power 
spectrum in that direction, potentially providing additional in¬ 
formation. Cross-correlation could also be useful in this re¬ 
gard, whether with intensity maps of other lines or galaxy sur¬ 
veys. The latter would include only the brightest galaxies but 
nevertheless trace the same cosmological structure. 

It is worth noting that, although we have focused this anal¬ 
ysis on CO, this same approach could just as well be applied 
to, say, [Cll] intensity mapping. [Cll] is also thought to trace 
star formation in high-redshift galaxies, but through a differ¬ 
ent (ionized) phase of the ISM. 

4.4. Astrophysical Contaminants and Complexities 

In this paper, we have assumed for simplicity that the cos¬ 
mological CO(l-O) map been perfectly isolated from the raw 
signal. Reality is certainly more complex. 

Most of the astrophysical foregrounds and backgrounds are 
expected to be continuum emission. These include the CMB 
itself, synchrotron radiation within the Milky Way, and (red- 
shifted) thermal dust emission from all galaxies along the 
line of sight. These components are expected to be relatively 
smooth in frequency space, so if the instrument has enough 
spectral resolution, the continua can plausibly be subtracted 
from the signal, while the remaining fluctuations should con¬ 
tain the cosmological CO signal. Note that this continuum 
subtraction effectively subtracts the lowest-k modes (k = 0 is 
the mean intensity), and so information in those modes will be 
reduced or removed. If those particular modes are important 
for the particular analysis at hand, the survey would need to 
be expanded to a larger volume to recover them. 

Line broadening and redshift space distortions, while not 
contaminants, will alter the power spectrum along the line 
of sight. This and previous studies have assumed a delta- 
function line profile, with zero width. However, if CO is a 
reasonable dynamical tracer within galaxies, then line broad¬ 
ening is an important consideration. We would expect broad¬ 
ening to reduce small-scale (high-k) power along the line of 
sight, limiting the narrowest useful frequency channel. For 
reference, a rest-frame width of 100 km s _l for CO(l-O) from 
Z ~ 2.5 is an observed frequency width of 11 MHz (roughly 
the frequency resolution of the “full” experiment in this paper, 
compare 8.25 MHz ~ 1 Mpc from §2.4). Additionally, galaxy 
peculiar velocities will cause redshift space distortions, alter¬ 
ing the apparent clustering of CO emission. Properly account¬ 
ing for these effects requires the cylindrical power spectrum 
P(k± ,ky) instead of the spherical P(k) in this study. 


The signal may be partially contaminated by interloping 
non-CO spectral lines that have been redshifted into the ob¬ 
served frequency band. While we expect low-/ CO lines to 
be the brightest lines in their spectral neighborhood, other 
lines may cumulatively add spurious fluctuations on top of 
the CO signal. In particular, HCN (z/ rest = 88.63 GHz from 
z ~ 1.6 —2.0) may contribute non-negligible shot noise to the 
measured power spectrum. In a 2D Q analysis, Breysse et al. 
(2015) found that using pixel masking to mitigate HCN fore¬ 
grounds resulted in a loss of shot noise information. As infer¬ 
ences about galaxy populations are sensitive to both clustering 
and shot noise components, this could complicate the astro- 
physical interpretation of an intensity mapping measurement. 
We note that, while data for the intensities of contaminant 
lines are currently sparse, in ALMA observations of lensed 
high-redshift star-forming galaxies, HCN is dimmer than CO 
by at least an order of magnitude (Spilker et al. 2014). 

Our observed band should also catch CO(2-l) from z ~ 
5.8 —6.7. We do not expect it to dominate the CO(l-O) sig¬ 
nal. However, removing the CO(l-O) signal could give us a 
residual CO(2-l) signal, if this removal can be done robustly 
enough, e.g. by cross-correlating with galaxy surveys. This 
would allow one band to simultaneously study two redshift 
ranges. 

Note that we do not expect the bright 158 pm [Cll] line to 
be a contaminant. For the observing band considered in this 
study, an interloping [Cll] line would have to originate from 
z > 55, when the universe is < 50 Myr old. It is unlikely that 
[Cll] will be found in significant abundance at these redshifts 
(though it would be a remarkable discovery). 

Finally, one more concern is that since the CMB would 
have a higher temperature at high redshift, then as an observ¬ 
ing background it may reduce the apparent luminosity of the 
CO lines (Obreschkow et al. 2009a; da Cunha et al. 2013). 
The severity of this effect depends on the ISM properties of 
high-redshift galaxies: if the CO-luminous regions of the ISM 
are significantly warmer at high redshifts, the CMB is less of 
a concern, and there is some evidence that this may be true 
(Harris et al. 2012; Daddi et al. 2015). However, suppose 
the CO(l-O) luminosity is globally reduced by a factor 0.4, 
which is similar to or more pessimistic than the predictions 
in da Cunha et al. (2013) at z ~ 2-3. P(k) would be reduced 
by a factor 0.16, and recalculating the fiducial detection sig¬ 
nificance SNR tot (§3.1, Eq. 7) yields ~ 1.67 a and 49.6er for 
the pathfinder and full experiments, respectively. While these 
numbers are certainly smaller, they still suggest a signal de¬ 
tectable by the full experiment. Nevertheless, the effect of the 
CMB as an observing background remains a valid concern 
due to unresolved uncertainties about the high-redshift ISM. 

4.5. Prospects for Cross-correlation 

Because intensity mapping smoothes over emission from 
many galaxies, cross-correlation of the signal with galaxy sur¬ 
veys offers a promising method of isolating the signal from 
galaxies alone (e.g. Pullen et al. 2013) at a specific redshift. 
The underlying idea is not that CO would necessarily origi¬ 
nate only from observed galaxies, but rather that individual 
galaxies and CO would trace the same large-scale structure. 
As a specific example, the survey area of the pathfinder ex¬ 
periment is well matched to certain galaxy surveys, including 
the 2 deg 2 COSMOS field (Scoville et al. 2007), which con¬ 
tains extensive multiwavelength data. Cross-correlating the 
CO map with existing galaxies can help validate the signal 
and may also constrain the properties of undetected galaxies. 
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Within the same field, multiple intensity maps of differ¬ 
ent lines could be quite complementary. CO is particularly 
well suited for cross-correlation with itself since it emits in 
a known ladder of emission lines, so, for example, CO(l-O) 
emission observed at z/ obs should correlate with CO(2-l) emis¬ 
sion at 2/y ohs . In addition, cross-correlation between CO and 
[Cll] signals from the same redshift, as well as between CO 
and Hi maps, can provide crucial information on the nature of 
various phases of the ISM in galaxies. 

Simultaneous intensity maps of Hi and CO emission, in 
particular, would be especially useful in probing the reion¬ 
ization epoch (e.g. Lidz et al. 2011). The former would trace 
the process of reionization, while the latter might trace star¬ 
forming galaxies. If reionization is in fact driven by star¬ 
forming galaxies, cross-correlating the two maps could place 
interesting constraints on the morphology of the reionized 
intergalactic medium (via Hi), the nature of the reionizing 
sources (via CO), and the relationship between the two. 

5. SUMMARY 

We have performed a preliminary analysis of the ability of 
CO intensity mapping to probe high-redshift galaxy popula¬ 
tions. The following list summarizes our main results: 

1. Based on our fiducial assumptions, we find that 
the CO(l-O) signal from z ~ 2 4-2.8 should be de¬ 
tectable by a realistic intensity mapping experiment 

(§3.1). Our model is consistent with a range of previ¬ 
ous predictions, but all models carry significant uncer¬ 
tainties, which are not definitively resolved by current 
observations. 

2. Details of the L co -M h relation produce measurable 
differences in shape of the power spectrum (§3.2). 
This encodes information about the galaxy population 
in the relative strengths of clustering (low-0 and shot 
noise (high-/:) modes. This implies that sensitivity to 
fluctuations over a large dynamic range of scales pro¬ 
vides better constraints on the underlying galaxy popu¬ 
lation. 
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3. We demonstrate the extent to which an intensity 
mapping observation can constrain the properties of 
galaxy populations (§3.3), namely (1) the L IR -L C o re¬ 
lation and (2) the CO luminosity function. This has sig¬ 
nificance for understanding molecular gas and its con¬ 
nection to star formation in high-redshift galaxies, to 
the extent that Lir traces SFR and Leo traces molecular 
gas. At high redshifts, those connections will require 
further study. 

We reiterate that intensity mapping directly probes the faint 
end of the luminosity function while similar inferences from 
galaxy surveys rely on extrapolation from observations of 
bright galaxies. This suggests that intensity mapping and 
galaxy surveys could be complementary avenues for under¬ 
standing galaxy populations at high redshifts. 

In summary, we have shown that intensity mapping is a 
promising method of observing CO in high-redshift galaxies. 
If CO traces molecular gas in these galaxies, then CO inten¬ 
sity mapping is a potentially informative probe of star forma¬ 
tion and molecular gas in the high-redshift universe, particu¬ 
larly in galaxies that will be hard to observe individually for 
the forseeable future. 
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Fig. 11.— Mean CO brightness temperature as a function of minimum halo mass, at z = 2.4. As M mln increases. (7co) naturally decreases, since fewer halos 
contribute to total CO luminosity. The shaded gray range indicates the values spanned by the (7co) spanned by the itlcr posteriors on SFR(M, z) from Behroozi 
et al. (2013a). The values of (7co) at /V/ nml = 10 9 , 10 10 , 10 11 , and 10 12 Mg have been labeled on the right. Note that (7co) was analytically calculated from the 
mass function of Sheth et al. (2001). 
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APPENDIX 

A. CO LUMINOSITY MODEL: SANITY CHECKS 
A.l. Minimum CO-luminous Halo Mass 

Here, we provide some justification for minimum CO-luminous halo mass of Mco, min = 10 10 M 0 in our model. From a practical 
standpoint, this informs our choice of dark matter simulation, which is not complete below ~ 10 1 h Mq. 

Physically, the cutoff is motivated by the idea that smaller galaxies at higher redshifts may not be CO-luminous. These galaxies 
may be less chemically evolved and/or less dusty than their low-redshift counterparts. As a result, CO in these systems might be 
less abundant and/or more likely to be dissociated by ionizing radiation. Such galaxies might not have significant CO emission 
despite active star formation. 

A strict mass cutoff is chosen for simplicity and is unlikely to be the most realistic model. However, the effect of loosening it 
is well within other model uncertainties, so the cutoff we have chosen is acceptable for the scope of this paper. 

To illustrate this quantitatively. Figure 11 shows the mean CO brightness temperature, (Tco), as a function of this cut¬ 
off mass Mco, min- Here, (Tco) has been obtained from the mean CO luminosity per volume, calculated as dLco/dV = 
Lco(M)(dn/dM)dM and converting to brightness temperature, where dn/dM is the analytic form of the halo mass 
function from Sheth et al. (2001). In our fiducial model at 7 = 2.4, lowering Mco. min does increase the mean CO brightness 
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Fig. 12.— CO line luminosity vs. galaxy stellar mass. Stellar masses were calculated from the stellar mass-halo mass relation of Behroozi et al. (2013a), the 
same source for our SFR-halo mass relation. For comparison, the observations of Genzel et al. (2010) and Tacconi et al. (2013) have been plotted: squares are 
measured values, triangles are upper limits. Note that the observed Leo data, in their respective studies, were converted to CO(l-O) luminosities from 3-2 and 2-1 
lines, assuming a certain scaling between the different line luminosities. There is significant uncertainty in this scaling. 

temperature, but with diminishing returns: Lqo(M) decreases more than the abundance of low-mass halos increases. To wit: in 
decreasing Mqo. m in by 1 dex from 10 11 M 0 to 10 10 M 0 , (7co) increases by ^25%. In further decreasing Me o, min from 10 10 M 0 
to 10 9 M 0 , (7co) decreases by less than ~ 10% — well within the uncertainty from SFR(M,z) alone — at M mm = 10 10 M 0 . 

Note that this particular justification is less valid if we assume Leo oc M, as in some previous models (§2.1 or Table 3). 
Incidentally, the fiducial models of those studies all assume Mco.min < 1O 9 M 0 , allowing more of the halo population to be CO- 
luminous. Note that these models do not take into account the rapid fall-off of the star formation efficiency below M* that is 
indicated by several empirical studies (e.g. Behroozi et al. 2013a and references therein). For comparison, the same models from 
Figure 2 have also been plotted in Figure 11 . 

This justification is also more uncertain at higher redshifts. There, for a fixed halo mass, the halo mass function grows steeper 
with redshift, which may balance out the decrease in LeoiM). 

A.2. Comparison with Existing Leo-Mass Observations 

Though sparse, measurements of Leo exist for galaxies at these redshifts. Comparing the modeled Lco(M) relation directly with 
observations is difficult, mainly because masses of dark matter halos are difficult to measure accurately. However, measurements 
of stellar mass are more commonly available from photometric data, and this is something we can compare. 

Figure 12 shows CO luminosity vs. stellar mass, both as seen in observed galaxies and as predicted by our model. Observa¬ 
tional data shown are from Tacconi et al. (2013) and Genzel et al. (2010) for z ~ 1-3 galaxies. Our model does not explicitly 
predict stellar mass, but we have calculated it using the stellar mass-halo mass relation from Behroozi et al. (2013a), which 
self-consistently yielded the SFR-halo mass relation used in our modeling. 

The general consistency between model and data is reassuring, though not surprising, since the Lir-Lco relation that went into 
the model was empirically fit to existing observations. What bears the most emphasis is the near-total lack of data points for 
galaxies with low stellar mass (< 1O 1O M 0 ), which, by number, comprise the overwhelming majority of the galaxy population. 
Note also that our model for the average properties includes quenched galaxies, while the data include only those galaxies detected 
in CO. 


B. IMPLEMENTATION DETAILS 

For clarity, this section describes additional details for calculating CO brightness temperature and the power spectrum. 


B.l. Halo Luminosities to 3D Intensity Map 

Once Leo is calculated for all halos, they are binned on a grid according to their positions. This grid may be defined in 
observational units (RA-Dec-frequency) or comoving units. Converting between them requires us to specify a cosmology and 
rest-frame line frequency. In the two dimensions across the sky, we set the grid resolution to be 10 times finer than the beam size, 
i.e. the voxel length in the RA and Dec directions was $fwhm/10. 

On a comoving grid, the observed brightness temperature from a voxel of comoving volume V vox , containing total luminosity 
Leo, vox, emitted at rest frequency ;/ rcst , from redshift z is 


c 3 (l+z) 2 Leo ,vox 
8nk B iy^ st H(z) Kox 


( ^rest 

GHz 


H(z ) 



= 3.1 x 10 4 /iK(l +z) 2 


km s -1 Mpc -1 


-l 
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B.2. Fourier convention and power spectrum 

We calculate the power spectrum through the Fourier transform, F(k), of the 3D CO temperature cube. The power spectrum 
P( k) (specifically, power spectral density) is then 


p(k) = v; u ' n \F(k )\ 2 


(Bl) 


where V surv is the comoving volume of the survey. 

The Fourier transform requires a choice of convention, which does not matter as long as it is consistently applied. However, 
because conventions and numerical implementations may differ slightly, we state ours here. 

For the continuous Fourier transform F( k), we use the following (non-unitary angular frequency) convention: 

Forward: F(k) = / f(x)e~ ikx d 3 x 


with the reverse transform being /(x) = (2w) 3 f F(k)e' k ' d 'k. We only need the forward transform for P(k), so this choice 
conveniently avoids additional factors of 2n. Our discrete Fourier transform F k is defined as 


N -1 

Fk = ^2 a„ exp 

n 



k = Q,...,n—l 


(B2) 


and is thus related to the continuous Fourier transform F(k) as F(k) « AxAyAz /'k, where Ax, Ay, and Az are the comoving 
voxel lengths. 


B.3. MCMC Implementation 

We used the emcee Python package (v2.1.0, Foreman-Mackey et al. 2013), an implementation of an affine-invariant ensemble 
algorithm (Goodman & Weare 2010), to perform the MCMC analysis. This approach prioritizes ease of implementation over 
computational speed, which we believe to be justified for the scope of this work. To sample the five-dimensional parameter space, 
we used an ensemble of 200 walkers taking 2500 steps (500,000 samples in total), after a fairly generous burn-in phase of 500 
steps. 

At each sampled position in parameter space, we calculate a power spectrum P mo d e \(k) by generating a new intensity map from 
one of 100 realizations of the survey volume. Because of the finite volume of these realizations, P mo d e \ is not immune to sample 
variance, so it carries its own error bars, a> model . 

The "observed" intensity map, in reality, would appear smoothed by the telescope beam (see Figure 1), i.e. convolved with 
a 2D Gaussian filter in each frequency channel. This attenuates the apparent power spectrum, especially at scales smaller than 
the beam. In practice, the MCMC procedure would require this convolution to be repeatedly computed, which is moderately 
expensive even with Fourier transforms. Instead, we calculate the "unsmoothed" power spectrum for P obs and P mol ic\ but increase 
the error bars on both by a factor [W(k)] _1 (see Eq. C4). Thus, the resolution limits are encoded in the error bars (crp model and 
o>„ bs ), rather than in P mode i and P obs . 

Within additive constants, the log-likelihood expression is 



[Pmode\(k)-Ppbs(k)f 

op^ao+oj.^m 


+ln [cr 


2 

^model 


(*)+4 obs 



(B3) 


where cr Pobs is given by Eq. Cl, and crp modd = B mo d e i/VMnodes /W, where N mo d es and W are functions of k given by Eqs. C3 and C4, 
respectively. 

The 3D Fourier transform makes calculating the power spectrum the most computationally expensive part of each MCMC step, 
particularly as dynamic range increases (larger survey volumes and finer resolution). 

C. ERROR BARS ON THE POWER SPECTRUM 

In this section, we summarize the calculation of o>(k), the uncertainty in the spherically averaged power spectrum, or the error 
bars on P(k )—and by extension, A 2 (k). This is essential for quantifying how precisely an intensity mapping experiment can 
detect a signal, as well as cosmological structure within that signal. 

The details that follow are previously mentioned in, e.g., Lidz et al. (2011) and Gong et al. (2012). We find that a small but 
non-negligible correction in the calculation is necessary (noted below). Otherwise, we have mainly redescribed the procedure 
here for completeness and clarity. 

This calculation accounts for three sources of uncertainty in the power spectrum: 


1. Sample variance. 

2. Thermal noise variance, from the instrument. 


3. Limited resolution. 
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In a real experiment, there can certainly be additional systematics to consider, but this calculation is an appropriate starting point. 
The equation that summarizes the calculation is 


crp(k) = 

Sample variance Thermal noise Resolution 

where 

• P(k) is the (observed) spherically averaged power spectrum, as a function of k 

• P n (k ) is the thermal noise power spectrum: 

P a (k) = cj 2 n V vm (C2) 

• N mo des (k) is the number of measured modes at k in a bin of width A k: 

Anodes (k) — ^ ^ (C3) 

• W(k ) is a factor that attenuates the power spectrum, arising from limited spatial resolution and mainly affecting high k beyond 
resolvable scales: 

W(k) = dp (C4) 

Jo 

We explain Eqs. C2, C3, and C4 in the subsections below. Before we do, the following points may help to clarify notation, 

language, and assumptions: 

• Redshift is z, at which H(z ) is the Hubble parameter, R(z) is the comoving radial distance, and D/iz) is the luminosity distance. 
We assume a single redshift 1 +z = t/ rest /I'obs = A 0 bs/A rest for the volume, an acceptable approximation for the surveys being 
considered. 

• Comoving Cartesian coordinates are labeled (x t . xi, X 3 ) in order to avoid confusion with redshift. 

• Xj and X 2 are perpendicular to the line of sight (denoted _L, mapping to RA/Dec). x 3 is parallel to the line of sight (denoted j|, 
mapping to redshift or observed frequency). 

• The instrument has system temperature T sys and number of observing feeds /V fcc ds • The full bandwidth of the spectrometer is 
A;/, centered on an observing frequency ;/ ( ,b s and divided into frequency channels of width 5 V . The observing beam has a 
Gaussian profile, with full-width half-maximum angle (?fwhm (obeam = 0FWHM/V81n2). The total observing time on the survey 
area is T obs . 


m 


Pn(k) \ 1 


VN modes (k) VN modes (k) J W(k ) 


(Cl) 


• The smallest resolvable 3D volume element is approximated as a cubic “voxel,” while the sky projection of all voxels along a 
line of sight is a “pixel” (subscripted “vox” and “pix”). V vox is the comoving volume of a voxel, while fl p ; x is the solid angle 
on the sky subtended by a pixel. The expressions for each are 


O — fT~ 

^ £pix — °beam 


and 


^vox [fl(z)ob earn] 


c 5„(l+z) 2 
H(Z) V rest 


(C5) 


• Vsurv ~ T 1 T 2 C 3 is the comoving volume of the survey, while O surv is the solid-angle sky coverage of the survey. 


C. 1. Noise Power Spectrum [P n (k), Eq. C2 ] 

Consider any voxel in the observed 3D temperature cube, covering a sky area fl p j x and frequency bin Sv as previously described. 
We assume that the random thermal noise fluctuations, i.e. a white Gaussian noise spectrum. By the radiometer equation, the rms 
temperature fluctuation per voxel is 


cr„ 


1 sys 


\jN feeds ^ T pi 


(C6) 


where r p j x = r 0 b s (fl p i x /fl surv ) is the observing time per sky pixel. 
The averaged power spectrum of this noise is a constant: 


Pn(k) = a„ 2 w 
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C.2. Number of Modes [N modes (k), Eq. C3] 

Because P(k ) is obtained from averaging P{ k) over a finite number of modes (points in fc-space), we need to know the number of 
independently measured modes used to compute the average. We label this number A mo d es and thus a factor 1 / \/N nmdcs ultimately 
enters into the calculation of ap(k ). 

Written explicitly. Abodes depends on both the scale of the fluctuation, k, and the choice of bin width, A k. Consider a radial 
k-space bin (i.e. spherical shell) between k and k + Ak. Assuming A k is chosen, the number of modes within this shell is 

N modes (k) = n(k) ■ 4nk 2 Ak (Cl) 


where n(k) is the number density of modes in k-space. Recalling the similar “density of states” calculation in statistical mechanics, 
n(k) is 


n(k) = 


/ 2n 2it 2n 

\U~l 2 T 3 


-1 


Esurv 

87T 3 


(C8) 


However, the power spectrum of a real-valued function is symmetric about k = 0, so only half of the modes contain independent 
information. This introduces a factor 1 /2 into the final expression, which is (from combining Eq. C7 and C8): 


Nmodcf k- A k) 


k 2 AkV SUTY 
4n 2 


C.3. Resolution Limits [W(k), Eq. C4] 

Finite spatial resolution means that information from high-k modes is lost, due to smoothing of features smaller than the beam 
or channel width. The basic reasoning here is that if the smoothed power spectrum P sm (k) results from attenuating the “intrinsic” 
power spectrum P(k) by a factor W(k) after smoothing, the uncertainty should scale as op(k) oc P(k)/P sm (k) = 1 /W(k). 

With arbitrarily fine resolution, the temperature field is T(xi,x 2 ,xf), which has a power spectrum P(k). In reality, we measure 
a smoothed field, T sm (x\ ,X 2 -X 3 ), which is the original field convolved with a Gaussian in each direction: 

T sm (x) = T(xi,x 2 ,x i )*G(xi\a l )*G(x 2 \a 2 )*G(x 3 \a 3 ) (C9) 


where G(x i |(Ti) = (27rCTi) _1,/2 exp[-xf/(2iT^)] (similarly for x 2 andx 3 ) is a normalized Gaussian function, and indicates convo¬ 
lution. In reality, the beam profile may not be perfectly Gaussian, and the frequency channels would be discrete bins. However, 
the approximation will suffice here. 

Angular resolution determines 0 \ and cr 2 , while frequency resolution determines o 3 . Specifically, 


t7 1 — tJ 2 — (T_j_ — R(z) CTbeam 


and 


& 3 = CTj| 


c S„(l+z) 

H (z.) t'obs 


(CIO) 


Since a convolution in real space is simple multiplication in Fourier space, the Fourier transform and power spectrum are 


T sm (k) = T(k u k 2 ,k 3 )e 




and 


P sm (k) = P(h , k 2 , k 3 )e-W£ + ^£ + fy® 


(Cll) 


where G = = 1 /cr j_ and Ti = 1 /a \\. The Fourier convention used here ensures the transformed Gaussian carries no normalization 

constant. 

We note a small but consequential correction: Lidz et al. (2011) and some following studies state 1; j = 27r/[R(z)<Tbeam]. but we 
find from the calculation above that it should actually be £j_ = 1 /erj_ = 1 /[R(z)(Jbeam] (without the factor of 27r). 

To get the spherically averaged power spectrum, define the variable // = cos i), where 0 is the spherical polar angle in k-space. 
Then we define components of k parallel and perpendicular to the z direction, k and k±, as 

ky = A '3 = kn and k± = sfk 2 -k 2 \ = k\/l- /f 1 (C12) 

and the expression for the power spectrum becomes 

P sm (k) = P(k ± ,k ll )e~^ + ^ ) 

= R(k,/r)^ 2 ^e Vi2 ( CT "“ CTi ). (C13) 


To get the spherically averaged power spectrum, this needs to averaged over 0 < /t < 1 (the upper half-space, since only half the 
modes are independent) at fixed k. In this paper, we assume P(k. //) is isotropic, so only the exponential needs to be averaged: 


P s Jk) = P(k) e 


-k 2 a 2 1 


e k V' 7 " a± )^ d[i. 


W(k) 


(C14) 


The expression after P(k) is W(k), which can be integrated numerically. 
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D. OPTIMIZING OBSERVING PARAMETERS 


The expression for ap(k) allows us to draw some conclusions about the effect of varying survey and instrument parameters. 
See Breysse et al. (2014) for a similar analysis with Q coefficients, i.e. a 2D power spectrum. 

We ignore the factor W(k) for this analysis, since its main effect is to limit the maximum k (smallest scale) at which at which 
meaningful measurements can be made. Certainly, it will also slightly suppress power at low k (large scales), but this should 
not be a dramatic effect. Then, the expression for crp(k), showing explicit dependences on instrument and survey parameters, is 
ultimately 


CTp(k) = 


2n 

k\J A k Ao 




Ul\Z) Jy feeds Aabs 


Test-“surv 




- D,.(z) 

'vm 


(Dl) 


where Di(z) = (1 +z)R(z) is the luminosity distance. Note that redshift z is directly determined by u 0 y, s . 

Written in this form, it is easier to see how instrument and survey parameters ultimately affect sensitivity. At a given scale k , 
the following adjustments will decrease a P (k ): 


• Increasing bandwidth: ap oc 1 / \/ Au. This expands the survey volume, increasing the density of modes in the k direction 
and thus Amodes- However, the observed redshift range will eventually grow large enough that cosmic evolution across the 
volume becomes important. 

• Decrease instrument noise: if ap is noise-dominated, ap oc T s 2 s /AfeedsTbbs- In that case, decreasing 7' sys , increasing Af ee ds, or 
increasing r Q b s will all serve to increase sensitivity. Eventually, though, decreasing system noise yields diminishing returns as 
sample variance becomes the dominant component of ap. 


• Optimize survey area. For a given mode k, the optimal survey area results in equal sample variance (oc fl surv 1 / 2 ) and noise 
variance (oc Q SU rv 2 )- In steradians, this is 

1 Nf C cdsTobs H (z) 


a 


surv,opt — 


Arest A ys ^ L (Z )] 2 

From noise From 


P(k). 


(D2) 


Because O survopt oc P(k), the predicted optimal area is model-dependent. 

As noted in §3.1, we can maximize the total detection significance, SNR tot (Eq. 7), by optimizing the area. For the observing 
parameters in this study, we obtain 0.6 and 9.2 deg 2 for the pathfinder and full experiments, respectively. Certainly, these 
optimal areas differ from our adopted survey areas (2.5 and 6.25 deg 2 ), but the improvement in SNR tot is not dramatic (7.89 to 
9.48 for pathfinder; 144 to 146 for full). 

However, in practice, there are other factors to consider. The minimum survey area cannot be arbitrarily small, because the 
arrangement of multiple feed horns on the instrument images a pattern of multiple beams on the sky. This pattern must 
be scanned across the sky to adequately cover the survey area, and limitations on the scanning strategy effectively set the 
minimum size of the survey area. This motivated our estimate of 2.5 and 6.25 deg 2 as the survey areas for the pathfinder and 
full experiments, respectively. 

The maximum survey area can be as large as the entire sky, but increasing total area would reduce the amount of integration time 
per unit area. This can decrease sensitivity to high-/: (shot noise) modes, which contain important astrophysical information. 


It is worth noting that the following parameters do not, for the simple assumptions above, affect ap (neglecting secondary effects 
such as redshift evolution across the volume and the curvature of the boundary surface, which we expect to be subdominant for 
the volumes being considered): 

• Resolution. Assuming the spatial scales being probed are much larger than the resolution limit, refining the resolution does 
not decrease a P significantly. 


• Survey area aspect ratio. For a fixed rectangular survey area, its relative dimensions have no effect. As the length grows, the 
width shrinks. In /--space, the density of modes becomes sparser in one direction but greater in the other, and A mo des does not 
change. 


E. OBTAINING DATA FROM OBRESCHKOW ET AL. (2009) 

For the results of §3.3.1, we downloaded the lightcone data of Obreschkow et al. (2009b) from their publicly accessible 
database. We selected all galaxies within a central 6.25 deg 2 area, located at comoving distances between 5770 and 6285 Mpc, 
which encloses the redshift range 2.4 < z < 2.8 in our assumed cosmology. From each simulated galaxy, we obtained four 
quantities: RA, Dec, comoving distance, and CO(l-O) flux. On the query form page 12 , the exact SQL query given was 

12 http : //s-cubed . physics . ox . ac . uk/queries/new? 
sim=s3_sax 
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select ra, decl, cointflux_l, distance 
from galaxies_line 

where ra between -1.25 and 1.25 

and decl between -1.25 and 1.25 

and distance between 5770 and 6285 

Since CO(l-O) flux was provided as an obesrved velocity-integrated flux, S v (units of Jy km s _I ), we converted to an intrinsic 
brightness temperature luminosity, L r (units of K km s -1 pc 2 ) following the derived equation in the appendix of Obreschkow 


et al. (2009a): 



(El) 


From these galaxy luminosities, we generated CO intensity maps and power spectra as described in §2.2.5. 
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TABLE 3 

Previous CO Intensity Mapping Models 


Paper 

Redshifts 

CO luminosity model 

Righi et al. (2008) 

z ~ 0—10 

Calculating SFR: SFR assumed to be driven by ma- Converting SFR to Lco : Based on M82 observa- 
jor mergers. Stellar mass M* formed in a merger of tions (WeiB et al. 2005): 
total mass M = M\ +M 2 calculated as: 

Leo = 3.7 x 10 3 SFR 

* M\Mo 

M = 4 ri 1 = 3.5 x 10~ 2 

n m m m 

Merger rate was calculated from extended Press- 
Schechter formalism (Lacey & Cole 1993). 

Visbal & Loeb (2010) 

z > 3 

Calculating SFR: SFR proportional to halo mass. Converting SFR to Lco : Based on M82 observa- 
Mhaio- tions (WeiB et al. 2005): 

11 / 1 3/2 L co = 3.7 x 10 3 SFR 

SF* = 6.2xlOT n M halo 

Gong et al. (2011) 

z ~ 6 — 8 

Calculating Lco : Fit a relation of the form At z = 6, 7, and 8 respectively: 

_ L (M h , lc \ b / + M halo y rf L 0 = 4.3 xlO 6 , 6.2x10 s , 4.0 x 10 6 L 0 

^ M c ) \ M c ) b = 2.4,2.6,2.8 

to the results from semi-analytic modeling by M c = 3.5 X 10 11 ,3.0 X 10 11 ,2.0 X 10 11 Mg 

Obreschkow et al. (2009a). J = 28 3433 

Carilli (2011) 

z ~ 6—10 

Calculating SFR: Cosmic SFR density calculated as Converting SFR to Leo: Combining SFR-Lir (Ken- 
the value required to reionize the universe at a given nicutt 1998) and disk galaxy Lir-Lco (Daddi et al. 
redshift (Bunker et al. 2010; Madau et al. 1999). 2010) relations: 

Lco = 1-0 X 10 4 SL7? 

Lidzetal. (2011) 

z > 6 

Calculating SFR: SFR proportional to halo mass. Converting SFR to Lco: Combining SFR-Lir (Ken- 
Mhai 0 : nicutt 1998) and Lir-Lco (Wang et al. 2010) rela¬ 

tions: 

SFR = 1.7 X 10~ 10 M halo 

Lco = 3.2 x 10 4 SFR 

Note: as originally derived, Lco ex (SFR) 3 / 5 . Linear 
scaling adopted for simplicity. 

Pullen et al. (2013) 

NOTE. — For consistent cc 
for full parameterizations anc 

z ~ 2 — 4 

mparison, re 
details. See 

Calculating SFR: Model A: SFR proportional to Converting SFR to Lco: Same as Lidz et al. (201 1): 

(Mhalo) / ' L co = 3.2 x 10 4 (SFRf /s 

SFR -1.2x10 (Mhaio) Note: for Model A, the effective model is still 

Model B-. SFR calculated from empirically fit SFR Lco 3L halo . For Model B, power spectra were 
Schechter functions. (Smitetal. 2012) obtained by multiplying those of A by a factor 

(7co,b) 2 /(7co,a) 2 (~ 4.8 2 at z = 3). 

ations here were simplified using the authors’ given fiducial assumptions, but we refer the reader to the original papers 
figure 2 for a comparison of Lco(^haio) between some of these models. 













